% File: bestfit.m
% Author: Peter Horak (Spring 2012)
%
% Purpose: fit a linear combination of functions to x and y data points
% using the least-squares fit
%
% Input: function, x-values, y-values, plot_answer
%   FUNCTION -- a function handle for a function that takes a single number
%       (i.e. x) and returns a row vector representing a sum of functions
%       evaluated with this number (i.e. [1,x,x^2]).
%   X-VALUES -- a vector of the x-values of the data
%   Y-VALUES -- a vector of the corresponding y-values
%   PLOT_ANSWER -- optional argument, 1 means make a plot of the data and
%       model, anything else means don't make a plot (this argument
%       defaults to zero)
%
% Example input: bestfit(@(x) [1,x],[1,2,3,4,5],[2,2.5,3.1,3.5,3.9],1)
%
% Output: a row vector of the coefficients to go with the corresponding
% functions from the function argument that result in a model closest to
% the real data.
%
% Errors: the program does not check that the function provided throught
% the handle is valid. Also, for some reason the program only works with
% anonimous functions.

function [coefficients,Rsquared] = bestfit(funct,xs,ys,plotanswer)

% Check the number of arguments (and optional argument)
% Based on an example from the website:
% http://blogs.mathworks.com/loren/2009/05/05/nice-way-to-set-function-defaults/
if nargin > 4
    error('Usage: too many arguments.');
elseif nargin < 3
    error('Usage: too few arguments.');
elseif nargin == 3
    plotanswer = 0;
end

% Check that the x and y data vectors are valid
if max(size(size(xs)) ~= [1,2])
    error('Usage: xs cannot have more than 2 dimensions.')
end
if max(size(size(ys)) ~= [1,2])
    error('Usage: ys cannot have more than 2 dimensions.')
end
[rx,cx] = size(xs);
if cx ~= 1
    if rx ~= 1
        error('Usage: xs must be vector (row or column matrix).')
    else
        xs = xs'; % Transpose vector to make column vector if necessary
    end
end
[ry,cy] = size(ys);
if cy ~= 1
    if ry ~= 1
        error('Usage: ys must be vector (row or column matrix).')
    else
        ys = ys'; % Transpose vector to make column vector if necessary
    end
end
[rx,cx] = size(xs);
[ry,cy] = size(ys);
if rx ~= ry
    error('Usage: input vector dimensions do not match.')
end

% Find the least-squares solution 
for i = 1:rx
    A(i,:) = funct(xs(i));
end
equation = cat(2,A'*A,A'*ys);
solution = rref(equation);
[r,c] = size(solution);
coefficients = transpose(solution(:,c));

% Error analysis
% I reviewed the meaning of R^2 on the webpage:
% http://en.wikipedia.org/wiki/Coefficient_of_determination
SStot = zeros(1,rx);
SSerr = zeros(1,rx);
for j = 1:rx
    SStot(j) = ((ys(j)-mean(ys))^2);
    SSerr(j) = (ys(j)-sum(coefficients.*funct(xs(j))))^2;    
end
SStot = sum(SStot);
SSerr = sum(SSerr);
Rsquared = 1 - SSerr/SStot;

% Optional plotting
if plotanswer == 1
    xs2 = linspace(min(xs),max(xs),rx*10);
    ys2 = zeros(1,rx*10);
    for k = 1:(rx*10)
        ys2(k) = sum(coefficients.*funct(xs2(k)));
    end
    plot(xs,ys,'+',xs2,ys2)
    legend('data','model')
    title('least-squares approximation')
    xlabel(cat(2,'R-squared value = ', num2str(Rsquared)))
end